Runge-Kutta Methods for Linear Ordinary Differential Equations

Three new Runge-Kutta methods are presented for numerical integration of systems of linear inhomogeneous ordinary differential equations (ODEs) with constant coefficients. Such ODEs arise in the numerical solution of the partial differential equations governing linear wave phenomena. The restriction to linear ODEs with constant coefficients reduces the number of conditions which the coefficients of the Runge-Kutta method must satisfy. This freedom is used to develop methods which are more efficient than conventional Runge-Kutta methods. A fourth-order method is presented which uses only two memory locations per dependent variable, while the classical fourth-order Runge-Kutta method uses three. This method is an excellent choice for simulations of linear wave phenomena if memory is a primary concern. In addition, fifthand sixth-order methods are presented which require five and six stages, respectively, one fewer than their conventional counterparts, and are therefore more efficient. These methods are an excellent option for use with high-order spatial discretizations. Introduction We consider the numerical integration of large linear inhomogenous systems of ordinary differential equations in the form du = Au g(t) ( 1) dt where A is an M by M matrix whose elements depend on neither u nor t, and u and g(t) are vectors of length M. Such essentially autonomous systems arise in the numerical solution of partial differential equations (PDEs) governing linear wave phenomena after application of a spatial discretization such as a finite-difference, finite-volume, or finite-element method. Examples of such PDEs are the linearized Euler equations governing acoustic waves and the Maxwell equations governing electromagnetic waves. The elements of A depend on the PDE and the spatial discretization. The inhomogeneous term g(t) is associated with either a source term or the boundary conditions. In the context of wave propagation, the system of ODEs is oftenmildly stiff with theeigenvaluesof A typically lying nearthe imaginaryaxis. The systemof ODEs arising from the applicationof a spatial discretizationto a systemof PDEs can be very large, especially in three-dimensional simulations. Consequently, the constraintson themethodsusedfor integratingthesesystemsaresomewhatdifferent from those which have driven much of the developmentof numericalmethodsfor initial value problems. Due to their high accuracyand modestmemory requirements,explicit Runge-Kuttamethods havebecomepopular for simulationsof wavephenomena[5,6,7,15,17]. Thirdandfourth-order methodsrequiring only two memory locationsper dependentvariable are particularly useful [3,13,14]. This property is easily achievedby a third-orderRunge-Kuttamethod[14], but an additional stage is required for a fourth-order method [3]. Since the primary cost of the integration is in the evaluationof the derivative function, and each stagerequiresa function evaluation,the additional stagerepresentsa significantincreasein expense. For thesamereason, error checkingis generallynotperformedwhensolvingvery largesystemsof ODEsarising from thediscretizationof PDEs. Therehavebeenseveralattemptsto developefficientmethodsfor integratinglinear systems of ODEs[4,9,10,11]. Thebasicpremiseof thesemethodsis that themajorcostin evaluatingthe derivativefunction is in forming thematrix A and the vector g(t). In the application considered here, the simulation of linear wave phenomena, the matrix A is never explicitly formed or stored. Hence the methods previously proposed for linear systems are not appropriate for this application. It is well known that a Runge-Kutta method with p stages has an order of accuracy not exceeding p [1,2]. For p_<4, methods of order p can be derived with p stages. However, fifthand sixth-order methods require at least six and seven stages, respectively. Nine stages are required for seventh-order accuracy and eleven for eighth-order accuracy [1]. Since the cost for our application is roughly proportional to the number of stages, this represents a significant limitation of higher-order Runge-Kutta methods. Several authors have considered various approximations to reduce the number of stages and the storage requirements of high-order Runge-Kutta methods. Shanks [12] was able to develop schemes with a reduced number of stages by requiring only that the accuracy conditions be approximately satisfied. Zingg et al. [16,17] propose methods with low storage requirements which are of high order for linear homogeneous ODEs but second-order otherwise. A similar idea was proposed previously by Lorenz [8].